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Abstract 

The stress fields of dislocations predicted by classical elasticity are known to be unrealistically 
large approaching the dislocation core, due to the singular nature of the theory. While in many 
cases this is remedied with the approximation of an effective core radius, inside which ad hoc 
regularizations are implemented, such approximations lead to a compromise in the accuracy of the 
calculations. In this work, an anisotropic non-singular elastic representation of dislocation fields is 
developed to accurately represent the near-core stresses of dislocations in a-iron. The regularized 
stress field is enabled through the use of a non-singular Green’s tensor function of Helmholtz- 
type gradient anisotropic elasticity, which requires only a single characteristic length parameter in 
addition to the material’s elastic constants. Using a novel magnetic bond-order potential to model 
atomic interactions in iron, molecular statics calculations are performed, and an optimization 
procedure is developed to extract the required length parameter. Results show the method can 
accurately replicate the magnitude and decay of the near-core dislocation stresses even for atoms 
belonging to the core itself. Comparisons with the singular isotropic and anisotropic theories show 
the non-singular anisotropic theory leads to a substantially more accurate representation of the 
stresses of both screw and edge dislocations near the core, in some cases showing improvements in 
accuracy of up to an order of magnitude. The spatial extent of the region in which the singular 
and non-singular stress differ substantially is also discussed. The general procedure we describe 
may in principle be applied to accurately model the near-core dislocation stresses of any arbitrarily 
shaped dislocation in anisotropic cubic media. 

Keywords: dislocation, non-singular, anisotropy, atomic resolution. Green’s function, molecular statics 
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I. INTRODUCTION 


The response of dislocations to complex short and long-range internal stresses has long 
been known to govern a wide array of macroscopic material phenomena, including basic 
mechanical properties. Such stresses may arise from external loading, or by interactions 
with internal material defects, such as neighboring dislocations, point-defects, and other 
material interfaces. Investigations of the collective behavior of discrete dislocation ensem¬ 
bles by dislocation dynamics (DD) methods^”® have proven successful in modeling material 
deformation at the micro-scale. Due to the convenience of having closed-form expressions 
for the elastic fields, such simulations have primarily focused on isotropic materials. For 
anisotropic materials, the literature is sparse with such studies, due to the heavy computa¬ 
tional requirements of numerical integrations in the elastic fields'. While efforts have been 
made to increase computational efficiency®’®, computational times can still be more than two 
orders of magnitude longer for anisotropic versus isotropic simulations. Thus, such simula¬ 
tions are currently relegated to very short times and very few dislocations. This has limited 
large-scale simulations of dislocation ensembles to materials with low anisotropy. 

The most common measure of a cubic material’s anisotropy is through the Zener 
anisotropy ratio^®. A, which is dehned by a material’s elastic constants as A = C 44 /C', 
where c' = (cn — ci 2 )/ 2 . The physical interpretation is that of a ratio of the resistance to 
shear on {100} planes in (100) directions (C44) to the resistance to shear in (110) directions 
(c'). Materials having a Zener ratio of unity are taken as isotropic. Iron is a material 
well known to exhibit signihcant anisotropy (^4=2.43) which drastically increases as the 
a — 7 phase-transition temperature (912 °C) is approached^^. Due to its abundance and 
favorable mechanical properties, it still remains as the basis of today’s most technologically 
vital materials. Several recent studies have employed classical anisotropic elasticity theory 
and specifically illustrated the importance of anisotropy in discrete dislocation calculations. 
Using anisotropic elasticity theory, Dudarev et al.^^ showed that the relative values of their 
elastic free energy may account for the anomalous abundance of a(lOO) rather than a/2(lll) 
dislocation loops observed experimentally at high temperatures. Their analysis showed such 
a conclusion would not be reached via isotropic treatments. Through DD simulations, 
Aubrey et al.^® showed the effect of anisotropy on the equilibrium shapes of dislocation 
shear loops over a wide temperature range. Experimentally observed slip-system-dependent 
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loop features such as sharp corners were able to be reproduced, where not possible using the 
isotropic theory. Similar DD simulations of Frank-Read sources showed that the isotropic 
elastic approximation leads to large errors in their activation stresses, further increasing 
with temperature^^. Dependence on orientation was shown to play a very important role, 
specihcally in allowing certain conhgurations much lower critical stresses to become sources 
than isotropic theories allow, and more aligned with experimental observations. Thus, the 
modeling of dislocations in a-Fe cannot be accomplished via isotropic elasticity. 

While modeling via anisotropic theories has made recent progress, current implementa¬ 
tions, including those cited thus far, must still disregard material contributions near the 
core due to the singular nature of such theories. At the same time, it is well known that 
many important physical processes undergone by dislocations, including, among others, their 
intersection, pinning/de-pinning, and capture/emission of point-defects, can be accurately 
described only in terms of the processes occurring at the sub-nanometer scale adjacent to the 
dislocation core. To avoid the singularity issue, a core radius^^ must be prescribed around 
the dislocation, typically several Burgers’ vector magnitudes in radius or larger, within which 
processes may occur in a pre-dehned, ad hoc fashion, losing the atomic-scale resolution. In 
cases where heightened resolution on the atomic aspects near the core are required, such 
assumptions are not-acceptable, as a complete model of the complex interactions of dis¬ 
locations is possible only when all material contributions governing such phenomena are 
included. One example where the near-core held is crucial is in the study of creep in met¬ 
als. Creep is a phenomenon governed largely by the climb of dislocations in response to 
their absorption and emission of vacancies^®. This phenomena is further enhanced as the 
otherwise isotropic diffusion of vacancies in metals is made anisotropic in the stress helds 
of dislocations due to drift forces caused by gradients in their mutual interaction energy. 
This interaction energy may be expressed as the product of the dislocation’s stress tensor 
with the vacancy’s (constant) formation volume tensor, as Thus, to capture the 

absorption process accurately, the stress held near the core must be known at the atomic 
scale to a high degree of accuracy. If such resolution could be obtained, a lattice kinetic 
Monte Carlo simulation could be employed to accurately calculate point defect capture rates 
of arbitrarily shaped dislocations. Such accuracy is not possible via the singular theories, as 
the predicted stresses become too large and capture radii are required. 

We seek here to develop and employ a computationally-feasible continuum dislocation 
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theory that can account for full material anisotropy, and can provide accuracy down to the 
atomic scale. While the end goal of a reliable multiscale material model may be obtained 
by the coupling of continuum to atomistic methods^^, such approaches require a large com¬ 
putational overhead and cannot offer the degree of flexibility of a purely continuum-based 
approach. Recent progress in employing gradient elasticity^® of the Helmholtz type for the 
modeling of dislocation £elds^®“^^ has shown much promise for the heightened resolution of 
DD simulations in isotropic materials, as illustrated by Po et al.^^. The attractiveness of this 
approach is that it leads to a non-singular description of the stress and displacement fields of 
dislocations and requires only one additional material parameter. With the recent derivation 
of the anisotropic non-singular Green’s tensor function of gradient elasticity of Helmholtz 
type by Lazar and Po^®, in this work we explore its application to describing the near-core 
stress helds of dislocations in anisotropic materials, where singular theories are not valid. 
Choosing a-iron as a model anisotropic material, we perform molecular statics calculations 
of screw and edge dislocations in a-iron using a recently developed magnetic bond-order 
potential (BOP) by Mrovec et al.^^. An optimization scheme is then developed to extract 
the required length-scale parameter based on the obtained atomic stresses. While the focus 
here is on the development and verihcation of a non-singular anisotropic dislocation theory, 
knowledge of the acquired length parameter allows the non-singular Green’s function to be 
applied to the modeling of other elastic helds, such as those due to point defects. We also 
note in a more recent work by Lazar and Po^^, a method to treat materials with varying 
degrees of anisotropy beyond cubic was developed, where the Green’s function used in the 
current work is a special case of that more general theory. 

In what follows, hrst a review of the classical anisotropic continuum expressions for dislo¬ 
cation displacements and stresses is given in section H. We then show that these classical ex¬ 
pressions can be made non-singular by the direct substitution of the non-singular anisotropic 
Green’s tensor function. The equations used to calculate the atomic stresses from the BOP 
relaxations are also described. In section HI, a description of the atomistic simulations that 
were performed is given, including the methodology in which the non-singular characteristic 
length parameter, £, was extracted from such calculations. Results and Discussion of the 
atomistic and continuum calculations and of the spatial extent to which the non-singular 
theory is important are given in sections IV, followed by a summary and conclusions from 
the work in section V. 
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II. STRESS FIELDS OF DISLOCATIONS 


In this section we derive analytical expressions for the Cauchy stress generated by a 
dislocation loop in both classical and gradient anisotropic elasticity. These expressions are 
implemented numerically in section IV and compared to atomistic calculations. 

The classical elastic theory of discrete dislocations can be introduced in two equiva¬ 
lent frameworks. The first is framework of compatible elasticity, which was established by 
Volterra^® three decades prior to the proposed existence of crystal dislocations^^“^° and five 
decades before the proof of their existence by early transmission electron microscopy^^’^^. In 
his approach, Volterra studied a compatible elastic problem in which dislocations are repre¬ 
sented by constant translational discontinuities of the displacement field across an internal 
material surface. In this paper, however, we adopt the alternative framework of incompati¬ 
ble elasticity originally proposed by Kroner^^ and Mura^'^. The main kinematic assumption 
of this framework is that the displacement gradient Uij is split additively into an elastic 
distortion and a plastic distortion 

= A? + ■ (1) 


In this kinematic framework, dislocations are sources of plastic distortion (see Fig. 1), and 
a particular form of the tensor f3^ will be specified for both the classical and the gradient 
theory. Given a specific form of /3^, another tensor of fundamental importance in the theory 
is the dislocation density tensor ck, which is obtained as the negative curl of that is: 


aij — 




P 

im,k * 


( 2 ) 


A. Classical anisotropic elasticity 

In classical anisotropic elasticity the strain energy density of the medium, W, is expressed 
as a quadratic form of the elastic distortions, that is 

W = ^Cijkiefjsfi = (3) 

where is the symmetric part of and Cijki is the standard rank-4 tensor of elastic moduli. 
Given W and the aforementioned kinematic assumption, the equilibrium equation, in the 


6 


5 


K 



C = dS dij 


FIG. 1. The plastic distortion is concentrated on the surface S, which is bounded by the dislocation 
line C = dS. 


absence of body forces, becomes: 

\ 9^- / ~ ~ ~ Pkl )j ~ LikUk — Cijklf^kl,j ~ 0 (4) 

\ ij J j 

where aij = dW/defj is the Cauchy stress tensor, and Lik = Cijkididj is the Navier differential 
operator. Eq. (4) can be solved using the Green’s function method. For an inhnite medium, 
the displacement and stress fields read, respectively^^: 


Ui{x) = -CmnpqG^i^^^ix) * /3™(a;) (5a) 

~ Cijkl^lqsCmnpqGkm,n{^) * ^psi^) ■ 

In eq. (5), is the Green’s function of the Navier operator and the symbol * indicates 
convolution over three-dimensional space. The plastic eigendistortion corresponding to a 
Volterra dislocation extending over a surface S is: 

( 6 ) 

where 6 is the Dirac delta function, and bi (Burgers vector) is the displacement jump across 
S. Because b is constant, the dislocation density tensor turns out to be concentrated on the 
closed line C = dS bounding the surface S (i.e. the dislocation line): 

= -^jkml3fr^,k = x')bi dL) (7) 


7 









Finally, substituting eqs. (6) and (7) into (5), we obtain the classical expressions for the 
displacement and stress helds of a dislocation loop, respectively: 


Ui{x) — j CmnpqG^ 

(8a) 

^ ^ijkl^lqs^mnpqGX ^hp<\L^ . 

(8b) 


Notice that both these expressions involve the derivatives of the anisotropic Green ten¬ 
sor. The classical Green tensor for anisotropic materials was hrst obtained by Lifshitz and 
Rozenzweig^® and Synge^^ using the Fourier transform method, and it reads: 

where Lij{k) = Cijkikkki is the Navier operator in Fourier space, and n is a unit vector on 
the equatorial plane normal to r. The classical Green’s tensor function and its gradient 
are singular at the origin, and thus the singularity is inherited in all of the fundamental 
equations of classical dislocation theory. 


B. Gradient anisotropic elasticity of Helmholtz type 


Gradient elasticity of Helmholtz type^°’^^’^® is a simplihed version of Mindlin’s gradient 
elasticity^®’^®’^*^ with only one gradient parameter, where the strain energy density of the 
medium, W, is expressed as a quadratic form of the elastic distortion and its gradient: 

IV = ^Cijki/3,pii + (10) 


For details regarding the reduction of Mindlin’s general theory to this form, the reader is 
referred to the work of Lazar and Maugin^® and Lazar and In Eq. (10), the charac¬ 

teristic length, £, is the only gradient parameter to be determined by atomistic calculations. 
Owing to this particular form of the strain energy density, the equation of mechanical equi¬ 
librium, in the absence of body forces, reads: 



/ dW \ 
\d^ ) 



d 


(7 ij Tj 


ijm,m 


) j ^ijkldj Pkl} LikL'^k L^ijklP 


( 11 ) 

where Tijm = dV//defj,^ = i'^Cijkiefim = is the double stress tensor, and L = 1 — 

is the Helmholtz differential operator. Eq. (11) is formally similar to (4), with the 





difference that the composed Navier-Helmholtz operator LikL replaces the Navier operator, 
and L(3^i replaces Letting Gim be the Green’s tensor of the composed Navier-Helmholtz 
operator, this similarity allows us to easily hnd the solution of (11) as: 


Ui 


( 3 ?) Cjympq^im,n * ^mnpq^im,n * f^pq I ^mnpq^ ^ )^p (12a) 


Gjkl^lqs^mnpqC^km,n * LcXps Gjkl^lqs^mnpq^km,n * ®ps y Gjkl^lqs^mnpq^km,n{,^ ^ )^p 


JC 

(12b) 

Notice that in order to obtain the second equality in both (12a) and (12b) we have considered 
a plastic eigendistortion spread according to 


* G (13) 

so that = I3^° and La„ = aj,. 

In order to implement Eq. (12b) numerically, it is necessary to know the analytical ex¬ 
pression of the Green’s tensor function of the Helmholtz-Navier operator. This was recently 
derived by Lazar and Po^^ using the Fourier transform method and it consists of the convo¬ 
lution between the Green’s tensors of the individual operators. It reads: 

1 pi 

G.,(x)=G»,»G = ^j( L-/(k) e-«’/'d,d^, (14) 

where the constant 1/i term in front of the integral replaces the singular 1/r term in equa¬ 
tion (9). The gradient of the non-singular Green’s tensor function of the Helmholtz-Navier 
operator reads^^: 

-I p27r pi 

Gij,m{^) =L-/{k) (15) 

where q = cos 9 and dq = — sin6*d6* (see left side of hgure 2). In gradient anisotropic 
elasticity of Helmholtz type, both the Green’s tensor function and its gradient are non¬ 
singular. However, compared to their singular counterparts, equations (14) and (15) require 
an additional integral over the azimuthal angle in Fourier space. 


C. Atomistic calculations of stress 

In order to justify the comparison between continuum and atomic stresses, we shall now 
discuss the dehnition of the latter. The equivalent atomistic description of the continuum 
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Cauchy stress has for some time been the subject of debate'^^"'^^. Being a concept derived 
and understood at the continuum level, complications arise in the calculation of atomic 
stresses, primarily when hnite temperatures are involved. Such complications stem from 
uncertainties in the dehnition of the volumes enclosing the material point of interest and 
how to adequately average atomic forces and velocities both spatially and temporally. Using 
a statistical mechanics approach, where the stress is described by the local momentum flux 
through an enclosed surface, Lutsko^'^ developed a microscopic stress tensor suitable for 
molecular simulations. Cormier et al.'^®, expanded this approach and showed its validity 
near highly strained regions in the material by comparing with the solutions of classical 
anisotropic elasticity theory. While several additional schemes have also been developed 
to treat various calculations^®”^®, the virial stress, based on the virial theorem of gases by 
Clausius^®, remains the most common method of calculation. In this representation, the 
general per-atom stress tensor for an atom a, can be expressed as 

. N 


(a) J- 

CT-• =- 

n 




hm 


/3=1 


(16) 


where and are the mass, velocity, and position of atom a, and is the 

force vector on atom a due to atom f3. For bcc metals the atomic volume may be taken as 
= a®/2, and the summation is taken over N nearest neighbor atoms within a prescribed 
cutoff distance. For the static calculations (T= 0 K) performed in this work, the hrst term 
in (16) vanishes and thus no averaging schemes are required, as the equivalent Cauchy stress 
is directly obtained^®.The force vector in (16) is computed directly from the total 

energy, 17*°*, as 


where 




p(«)(/5) 


Qjjtot j.(a){l3) 
Qj‘(a){l3) ||]-( q :)(/3) 


(17) 


r(/5). We note that equation (16) is not momentum conserving and 


may produce erroneous results (e.g. non-zero normal stress components at free surfaces) 
when used in certain cases^^’'^’^. 


III. CALCULATIONS 

The anisotropic non-singular (A-NS) theory is enabled through the extraction of the 
characteristic length parameter which we obtain through a htting procedure with the results 
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of atomistic calculations. In this section we describe how these calculations were performed 
and the procedure by which the length parameter was subsequently obtained. 


A. Atomistic Calculations 

Molecular statics (MS) relaxations of 1/2(111) screw and 1/2(111){110} edge dislocations 
in iron were performed. The atomic interactions were treated with a recently developed 
magnetic bond-order potential (BOP) for iron by Mrovec et al.^^. Based on the tight- 
binding approximation, bond-order potentials'*^ offer a link between electronic structure 
and traditional molecular treatments. This potential also incorporates the Stoner model of 
itinerant magnetism to capture the effects of magnetic interactions between iron atoms. More 
commonly used interatomic potentials for iron such as that of Ackland et al.®^, struggle to 
reproduce important directional bonding characteristics and are not capable of addressing 
the issue of magnetism, which in iron has important consequences. Furthermore, such 
potentials have been specihcally known to poorly reproduce such properties as the Peierls 
barrier for the motion of 1/2(111) screw dislocations compared to density functional theory^^ 
and experimental results. In contrast, the present BOP has been shown to very accurately 
reproduce bonding characteristics and dislocation properties in iron. 

A schematic of the simulation cell used in the atomistic calculations is shown to the right 
in hgure 2 with orientation and geometric parameters listed in table I. The dislocations, both 
with line directions along the (periodic) 2 ; direction, were introduced in the center of the box 
by displacing all atoms according to the anisotropic singular (A-S) elastic solutions (equation 
(9)). Atoms in the hxed region were then held rigid, while atoms in the ’active’ region were 
allowed to relax according to the Fast Inertial Relaxation Engine (FIRE) minimization 
scheme^^. The relaxations were considered “complete” when the largest force on any atom 
is less than 0.01 eV/A. Stresses in this relaxed conhguration were then computed according 
to (16). 

B. Extraction of the Characteristic Length 

The characteristic length ^ in the gradient elasticity formulation equations (10), (14), 
and (15) may be obtained once the BOP relaxations have been performed. To accomplish 
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FIG. 2 . Schematic of the geometries used in the continuum stress calculations (left) and the 
atomistic simulations (right). 


X y z 

Lx/Cl 

Ly/a 

Lz/ Oj 

Active (Fixed) Atoms 

1 / 2 ( 111 ) screw dislocation [ 112 ] [ 110 ] [ 111 ] 

92 

92 

2 

11550(7344) 

1 / 2 ( 111 ){ 110 } edge dislocation [ 111 ] [ 110 ] [ 112 ] 

110 

108 

V6 

55056(3498) 


TABLE I. Orientation and geometric parameters used in the atomistic calculations (lattice param¬ 
eter a = 2.85A). 

this, a fitting procedure based on the minimization of an objective function is implemented. 
The objective function is the sum of the spectral norms of the residual matrices, denoted 
as SSNR(^\ over a chosen group of N (~16-24) atomic positions (r^")) surrounding the 
dislocation core. The residual matrix represents the difference of the calculated stress tensors 
from the non-singular continuum theory, and that from the atomistics, Thus, 

the objective function can be expressed as 

1 ^ 

SSNR(^) = (18) 

and is normalized by the shear modulus, fi, where /i = C44. The optimal characteristic length 
parameter can then be found as the value of i that minimizes (18). 

To better understand the role of each stress component in determining £, a matrix of the 
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sum of the norm of the residual components, denoted as was calculated as 







( 19 ) 


A plot of the components of SNRjj versus I will reveal the sensitivities and dependencies of 
I in reproducing the atomistic stresses, as will be seen in the section IV. 


C. Continuum Stress Calculations 

A schematic of the geometry used for the continuum calculations is given in the left side of 
figure 2. For comparison, calculations are also performed using the singular anisotropic (A- 
S) and isotropic (I-S) solutions from the classical approach. The calculations are performed 
at the reference atomic positions super-imposed onto a plane normal to the dislocation line. 
To calculate the stresses of an inhnitely long straight dislocation, we use the loop equations 
for stress (( 8 b) and (12b)), letting the loop radius go to inhnity. Thus we must choose an 
appropriately long dislocation line such that the contributions at and beyond the endpoints 
on the line are negligible. A dislocation line length of 306 was found to be sufficient, where 
6 is the magnitude of the Burgers vector. The plane where the calculations were performed 
was placed intersecting the midpoint of the dislocation line. The line was then populated 
with 50 evenly spaced points for the numerical integration of the stress equation. For each 
integration step along the line, integrations of the Green’s tensor function are required. 
These Fourier-space integrations are performed using Gaussian quadrature integration on 
the unit circle lying on the plane perpendicular to §3 (A-S) and over the unit hemisphere 
(A-NS). For the former, 8 quadrature points were found to be sufficient, while for the latter, 
the most critical calculation required 30 and 200 points for the integrations over 0 and 
g, respectively. The elastic constants obtained from the iron potential (cn = 243 GPa, 
C 12 = 145 GPa, and 044 = /! = 119 GPa) were used in all elasticity calculations. We note 
that for the isotropic-singular (I-S) solutions, the same methodology can be followed as the 
A-S with the substitution of Cn —)■ 2 C 44 -|- C 12 . 
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IV. RESULTS & DISCUSSION 


A. Screw Dislocation 


Figure 3 plots the results of the optimization of ^ for both screw (left) and edge (right) 
dislocations. Components of the sum of the normed residual (SNR^^^) matrix versus char¬ 
acteristic length parameter are plotted with inset plots showing the sum of the spectral 
norm of the residual matrices, SSNR*'^^ whose minimum determines the overall best choice 
of 1. For the screw dislocation, the minimums of the yz and xz components of the 
matrix give values of £ of 0.232 a and 0.316 a, respectively. In the inset plot of SSNR^^^ 
the minimum is found at a value of ^ =0.288 a, which we take as the best £t value for £. 
In figure 4, the results of the BOP relaxations and continuum calculations are plotted for 




FIG. 3. Plots of the components of the sum of the normed residuals (SNR^^^) versus characteristic 
length parameter for screw (left) and edge (right) dislocations. The minimum of each component 
gives the best value of (. considering only that component. The inset plots show the sum of the 
spectral norm of the residual matrices, SSNR^^\ whose minimum (encircled) determines the overall 
best choice of £, considering the contributions of all components. 

two rows of atoms below and two rows above the dislocation core in the [110] direction. In 
each plot, the abscissa represents the [112] direction in units of lattice parameters. It is 
immediately clear that the isotropic singular (I-S) solution gives a drastic over-estimation of 
the atomistic calculations, most notably in the xz component, where errors as high as 3-4 
times are observed. The A-S solution, however, yields stress values remarkably close to the 
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BOP calculations, within several percent, with the exception of the xz—component of the 
atoms belonging to the core. Results of the A-NS calculations are plotted as green points 
(where £=0.288 a) on a vertical line whose endpoints have horizontal bars representing the 
lower (black) and upper (red) values of £ = 0.232 a and 0.316 a, inferred from the range of 
minima of the components of SNR^y in hgure 3. The purpose of this line is to illustrate 
the sensitivity of the stress calculations to variations in the chosen £, which in this case we 
hnd to be negligible. For the screw dislocation, the A-S and A-NS solutions both provide 
excellent agreement with the BOP stresses for all atomic locations except for the core’s 3 
nearest neighbor atoms, where A-NS gives values of 0-22% larger than the atomistic calcu¬ 
lations and is thereby much closer to atomistics than the singular helds which deviate by 
65-103%. 


B. Edge Dislocation 

For the edge dislocation, the deviations between the components of the matrix 

in hgure 3 are found to be larger, ranging from a minimum of £ =0.386 a {zz component) 
to a maximum £ =0.860 a {xy component). While this represents quite a large range of 
candidate values for £, we hnd that due to its large sensitivity to changes in £ and its highly 
localized minimum near £ =0.4 a, the xx component of SNR*-^^ has the most inhuence on 
the minimum of SSNR*^^\ giving it a well-dehned minimum at the best ht value of 0.470 
a. The results of the BOP relaxations and continuum calculations are plotted in hgure 
5, for the two atomic rows below, and two atomic rows above the glide plane in the [110] 
direction. In each plot, the abscissa represents the [111] (Burgers) direction in units of 
lattice parameters. For the edge dislocation, the results show the BOP relaxations diher 
signihcantly from both of the I-S and A-S solutions. As expected, the I-S solutions again give 
a drastic over-estimation of each component of stress, in some instances nearly an order of 
magnitude too large. The A-S solutions oher some improvement, keeping the same trends of 
the I-S solutions, but decreasing their magnitude by typically between 25-50%. The trends 
predicted by the singular theories, however, are not observed in these atomic rows nearest to 
the glide plane. This is where we hnd the greatest beneht in implementing the non-singular 
solution. In hgure 5, the A-NS solutions are presented as green points (where £ = 0.470 a) on 
a vertical line that goes from £ = 0.386 a (black bar) to 0.860 a (red bar), taken from hgure 
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FIG. 4. Normalized stress calculations for a screw dislocation at atomic locations two rows above 
(y=\/2/4a, 3\/2/4a) and two rows below the dislocation core (Y=-\/2/4a, -3\/2/4a) in the [110] 
direction. In each plot, the abscissa, Y, represents the [Il2] direction, and the dislocation core is 
centered at X=0, y=0. Plotted are the singular anisotropic and isotropic elastic solutions (A-S 
and I-S), the results from the atomistic simulations (BOP), and the anisotropic non-singular elastic 
solutions (A-NS) using the best fit value oi £ = 0.288a. Horizontal black and red line segments 
represent A-NS stress calculations performed using the bounding values of £ taken from figure 3, 
as discussed in section IV A, and are connected with a vertical line to aid in visualizing. Black and 
red segments represent stress calculations using £ = 0.232 a and 0.316 a, respectively. 
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FIG. 5. Normalized stress calculations for an edge dislocation at atomic locations two rows above 
(y=\/2/4 a, 3\/2/4 a) and two rows below the glide plane (Y=-\/2/4 a, -3\/2/4 a). In each plot, 
the abscissa, X, represents the Burgers direction, and the dislocation core is centered at X=0, 
Y=0. Plotted are the singular anisotropic and isotropic elastic solutions (A-S and I-S), the results 
from the atomistic simulations (BOP), and the anisotropic non-singular elastic solutions (A-NS) 
using the the best fit value of £ = 0.470a. Connected horizontal black and red line segments 
represent A-NS stress calculations evaluated at -^ = 0.386 a and 0.860 a, respectively. 


3. At the two nearest rows to the glide plane, we see very good agreement with the BOP 
calculations, for every stress component. The very small-magnitude sinusoidal behavior of 
the xy component is also very nicely captured, where the singular solutions have failed. One 
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notable characteristic for the edge dislocation is the much larger variations in the calculated 
stresses over the range of candidate £’s used. This is most prominent in the hydrostatic 
stress components in the V = ±3\/2/4 a atomic rows. The A-NS solutions predict a weaker 
decay in stress than observed in the atomistics, suggesting improvement with a larger value 
of i. However, the differences compared to the BOP calculations are always relatively 
small compared to the singular cases which can be off by a factor of greater than an order 
of magnitude. A unique aspect worth noting was the ability of the A-NS calculations at 
multiple atomic locations to give the correct sense of the atomic stresses, where the singular 
solutions failed. This is seen in the atomic rows on the glide plane in ayy, and for some 
positions in axy 


C. Practical implications 


The results presented thus far have illustrated the drastic increase in accuracy in imple¬ 
menting the anisotropic non-singular solutions near the dislocation core. However, since this 
additional accuracy comes at a computational cost due to the additional nested integral in 
the non-singular Green’s function, in practice it is vital to know the spatial extent to which 
the A-NS solutions are necessary and where the singular solutions are sufficient. To do this, 
we have plotted in figure 6 a percent difference field of the A-S stress tensor with the A-NS 
stress tensor for a screw dislocation (left) and and edge dislocation (right), according to the 
expression 


%DIFF(r)^-s = 100% ■ 


kA-S(r)-aA-NS(, 

II^A-NSMII 


( 20 ) 


In the figure, an overlaid grid of points represent a reference atomic lattice for illustrative 
purposes. Additionally, all locations with a percent difference greater than 100 and 250 are 
lumped into the same color for the edge and screw dislocation, respectively. For the screw 
dislocation, the A-S and A-NS solutions converge to less than 20% and 3% difference at radial 
distances of approximately 0.9 b to 3 b from the core, respectively. For the edge dislocation, 
the differences extend to a much longer range in the material. Along the atomic rows adjacent 
to the glide plane, the A-S differs from the A-NS solution by > 100%, 50%, 20%, 10%, at 
|x| < 1.0 b, 1.9 6,3.3 b, and 6 b, respectively. Along a line perpendicular to the glide plane 
and intersecting the core, the A-S differs from the A-NS solution by > 100%, 50%, 20%, 10%, 
at \y\ < 0.65 b, 1.0 6,2.2 6, and 4.4 6, respectively. 
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FIG. 6. Percent difference field of the A-S stress tensor with the A-NS stress tensor (calculated 
from (20)) for a 1/2(111) screw dislocation (left) and a 1/2(111){110} edge dislocation (right). 
Corresponding reference atomic lattices are super-imposed for convenience. 

Similar calculations of the percent difference field of the I-S solutions with the A-NS 
show a convergence to a (roughly) constant difference of ~ 85% occurring just beyond 
the core atoms for the screw dislocation. For the edge dislocation, along the atomic rows 
adjacent to the glide plane, the I-S differs from the A-NS solution by > 100% and 50% 
at |a:| < 1.7 6,4.0 b, respectively. At larger distances, the differences oscillate between a 
minimum of ~ 20 % to a maximum of ~ 50 % for > 46. 


D. General discussion 

Implementation of the non-singular approach requires a priori knowledge of the char¬ 
acteristic length for the material. While in this work we have focused solely on near-core 
dislocation stress helds from atomistic calculations as means to extract this material pa¬ 
rameter, there are a variety of general approaches that one may follow to obtain it. Using 
density functional theory (DFT), Shodja et al.^^ were able to calculate strain-gradient char¬ 
acteristic lengths for several nearly isotropic materials {A 1), by relating the calculated 
Hessian matrix of the system to the double stress tensor. Other approaches include the de¬ 
velopment of analytical expressions derived from interatomic potentials'^. In our approach, 
we found that htting to the helds of both screw and edge dislocations did not result in a 
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single characteristic length. On the contrary, we found to be roughly 63% larger than 
j^screw^ However, noting the extremely minute variations of the stresses with changes in I 
seen in hgure 4, we can easily justify the use of the length parameter obtained from the edge 
dislocation (£=0.470 a) to describe both dislocations, and thus such a discrepancy poses no 
conflict in the case of iron. For other materials, this aspect remains an outstanding issue. 


Current methods of point defect and dislocation simulations stand to beneht greatly from 
the increased accuracy that the non-singular implementation provides. In DD simulations, 
ad hoc assumptions of the short range interaction of dislocations may now be replaced 
with non-singular helds based on and verihed by atomistic calculations. This includes a 
possible re-de£nition of the dislocation self force that may account for near-core contributions 
and material anisotropy, where the line-tension approximation is insufficient. Alternatively, 
simplihed descriptions of near core behavior in DD simulations may be calibrated to match 
the non-singular solution for the helds near the dislocation core. Mechanisms of dislocation 
junction formation as well as pinning and de-pinning and interactions with precipitates may 
also be better understood when such contributions are accounted for. In the same way as 
demonstrated in this work, the method could be implemented for the other transition metals 
and thereby help to analyze and possibly understand the differences of the bcc transition 
metals with respect to anomalous slip^®’^^. 


The diffusion of point defects in dislocation stress helds is known to be anisotropic due to 
drift forces acting on the point defect resulting from the interaction of the dislocation stress 
held and the point defect’s formation volume. In simulations of this type, the large overpre¬ 
diction of the dislocation stresses by the singular theories leads to large overpredictions of 
the interaction energy, and thus, important material parameters such as dislocation capture 
rates for interstitials and vacancies are inherently unreliable. Among other issues, this has 
important consequences for simulations of creep and swelling behavior in irradiated metals. 
In the case of the dislocation bias factor^®, a weaker interaction held may explain the order 
of magnitude overestimation of bias factors calculated with elasticity theory compared to 
those expected from empirical swelling data. 
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V. SUMMARY & CONCLUSIONS 


An anisotropic non-singular continuum theory of dislocation fields was derived and ap¬ 
plied to accurately represent the near-core stresses of dislocations in iron. The theory was 
enabled through the use of the recently developed non-singular Green’s tensor function of 
gradient anisotropic elasticity of Helmholtz-type^^, which requires only a single characteristic 
length parameter in addition to the material’s elastic constants. Using a magnetic bond- 
order potential to model iron interactions, molecular statics calculations were performed, 
and a htting procedure was developed to extract the optimum length parameter required by 
the non-singular theory. Results show the method can accurately replicate the magnitude 
and decay of the near-core dislocation stresses even for atoms belonging to the core it¬ 
self. Comparisons with the singular isotropic and anisotropic theories show the non-singular 
anisotropic theory leads to a substantially more accurate representation of the stresses of 
both screw and edge dislocations near the core. 

Conclusions of our study can be summarized as follows: 

• We have shown that an anisotropic non-singular expression for the stress as a line 
integral may be obtained by direct substitution of the anisotropic non-singular Green’s 
tensor function gradient of the Helmholtz-Navier type^^ into the classical line integral 
expression. 

• A htting procedure was developed and shown to be an adequate methodology for the 
extraction of the material length parameter, i, required by the non-singular theory. 
In bcc iron, this parameter is found to be £ = 0.288 a and 0.470 a when htted to the 
atomic stress helds of screw and edge dislocations independently, respectively. Due 
to the marginal sensitivity of calculated stresses on i for the screw dislocation, it is 
sufficient to take the value of 0.470 a, extracted from the edge dislocation calculations, 
as a single, best value for i. 

• For screw dislocations, the non-singular stresses match extremely well with the BOP 
calculations, even including the core atomic locations, where the singular solutions 
over-predict stresses by factors of greater than two. 

• For edge dislocations, the non-singular solutions offer a drastic improvement in match¬ 
ing the BOP calculations over the singular solutions, where stresses are found to be 
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as high as an order of magnitude too large in some cases. At multiple locations, the 
non-singular solutions are also found to match the correct sign of the atomic stresses, 
where the singular solutions give incorrect signs. 

• In practice, for an arbitrary dislocation loop or segment, it is advised that if less 
than 10%(50%) error with the non-singular solution is desired, the anisotropic-singular 
solutions should be implemented no farther than a cut-off radius of ~6&(1.96) from the 
core. For similar cut-off distances, the isotropic singular solutions oscillate between 
~20-85% error, depending on the angle relative to the core and the character (screw, 
edge, mixed) of the current segment being integrated. 

The results of this study highlight an advance in the capabilities of continuum elasticity 
approaches to resolve material stresses at sub-nanometer distances from dislocations, typi¬ 
cally reproducible only by atomistic simulations. The general procedure described may in 
principle be applied to model any cubic anisotropic material where high accuracy calcula¬ 
tions of near-core dislocation stresses are required. 
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